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§1. Introduction. Recently, Li and Sanders [2] have introduced a class of finite differ- 
ence schemes to approximate generally discontinuous solutions to hyperbolic systems of 
conservation laws. These equations have the form 

Q 

^q + V • g(q) + s(q) = 0 
q(0) = q 0 , 

together with relevant boundary conditions. When modelling hypersonic spacecraft reen- 
try, the differential equations above are frequently given by the compressible Euler equa- 
tions coupled with a nonequilibrium chemistry model. For these applications, steady state 
solutions are often sought. Many tens (to hundreds) of supercomputer hours can be de- 
voted to a single three space dimensional simulation. The primary difficulty is the inability 
to rapidly and reliably capture the steady state. In these notes, we demonstrate that a 
particular variant from the schemes presented in [2] can be combined with a particular 
multigrid approach to capture steady state solutions to the compressible Euler equations 
in one space dimension. We show that the rate of convergence to steady state coming from 
this multigrid implementation is vastly superior to the traditional approach of artificial 
time relaxation. Moreover, we demonstrate virtual grid independence. That is, the rate of 
convergence does not depend on the degree of spatial grid refinement. 

Before continuing, we review the particular variant of the numerical discretization of 
d*g(q) we wish to combine with multigrid. We assume that the problem to be solved is 
hyperbolic. That is, the Jacobian matrix of g(q), denoted here by d q g(q), has real eigen- 
values Aj(q), t = 1, . . . ,m sind a complete set of eigenvectors r;(q). The approximation of 
# z g(q) in grid cell j is given by 

^ (^*(qj+i/2;q;+i/2>qj+i/2) — kg(qj-i/2;qj-i/2>qi-i/2)) » 

where hg is a two point numerical flux function consistent to g. That is, hg(p; q, q) = g(q). 
The particular numerical flux function we take below is of Roe type [l] 

hg(p; a, b) = ^ (g(a) + g(b) - J2(p)|A(p)|X(p)(b - a)) , 

where J2(p) is the matrix of right eigenvectors to d q g(p), L(p) the inverse to R(p) and 
I a (p)I = v $)• We let q,+i /2 = §(qH-i + q,) and compute q} +1/2 and qJ +1/2 

according to the following recipe. 

(i) At each cell interface j + 1 /2, compute 

1)1,0 = i(<b+i/a)(q,’+i “ <lj) 

■0 2,+ = £(qj+i/ 2 )(qj +2 — 2q,+i + q,) 

U 2 -" = £(q i+ i /2 )(q j+1 - 2q,- + q^-i) 
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and componentwise (i — 1, . . • , to) 

_ /min(|I>J’"|,3|I>J’ 0 |)sigii(P?“) if Mqj+ 1 / 2 ) > 0 ^ 
(cj +1 / 2 )i - j min (|p?.+|,3|I>J’ 0 |)sign(D?’ + ) otherwise 

to obtain interface values 

qJ+ 1/2 = qj+i/2 - ^(q;+i/2) c i+i/2- 


(ii) Limit the interface values by 

<*5+1/2 = c b v 

+ J2(qj) ^minmod(L(qj)(qJ+i /2 — q^')*/*^ **/)(<*;— 1/2 <*i^v 

q i+i/2 = t L'+ 1 v 

+ R( qj+i) ^minmod(X(qj+i)(qj+i/2 _ qi+i)»P-^(fli+i)( t li+3/2 _ c b‘+ 1 ))j 
where p > 1 is called a compression factor. 

Note that for p taken larger than 1 we have generically that qj+ 1/2 — <*j+i/2 — c *j+ 1 /2‘ 
Therefore, generically 

^(qi+i/aJ^j'+i/a’^i+i^) = 8 ( qj+1/2 )• 

In the case when g(q) = q and q € R 1 , the finite difference formula above reduces 
generically to 

(Jjq)y « ^ (|(«u+i - %->) - j(^+> “ Sq J + 3q >-‘ “ qj “ 2> ) ‘ 

(See Section 3.2.) 

In the next section we motivate the basics of multigrid and study its convergence 
for an upwind scheme. In Section 3 a hybrid of multigrid with an approximate Newton s 
iteration is introduced and analysed. In the last section, we numerically treat ; ““JJSid 
supersonic flow problem using the difference formula given above together with a mul ign 

algorithm developed below. 
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§2.1 The basics of multigrid. Textbooks are filled with iterative techniques to solve 
large, sparse linear systems. Such systems are generally encountered when seeking ap- 
proximations to steady-state solutions of partial differential equations in one, two or three 
space dimensions. The rate of convergence offered by almost all iterative techniques un- 
fortunately depend on the size of the system. Specifically, by increasing the number of 
spatial grid points, the iterative technique will require a greater number of iterations to 
achieve a given error reduction. This phenomena can be exhibited by considering the finite 
difference scheme 


Ax 


(«;-«j-,) + /(*i) = o j = 1,2,..., J, 


Ax 


= j, 


with periodic boundary conditions, 

u 0 = 

used to approximate the solution to the differential equation 

+ ,{z) = °’ 
u(0) = u(l). 

(The finite difference scheme above has a nonunique solution provided that f(xj) = 

0.) Letting u represent the J dimensional vector of unknowns Uj, and f the J dimensi onal 
vector of /(ij), we can write the finite difference scheme above symbolically as 


( 2 . 1 ) 


Dju + f = 0, 


where Dj is a J x J matrix, and we will consider the rate of convergence of the artificial 
time relaxation technique 


( 2 . 2 ) 


u fc _ u fc-l _ + 


U° = V, 


where p > 0 is a relaxation parameter and v is an arbitrary starting vector. Clearly, if 
lim*_oo u* exists, the limit solves (2.1). 

Recall the discrete Fourier transform of a J-periodic grid function vj: 




and its associated inversion formula, 
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ut u* denote the kth iterate coming tom (2.2), let u denote the unique stationary solution 

to (2.1) normalized so that ZU “1 = “?■ “ d “ “ d '” 0 ‘' ‘ h ' ^ 

time iteration error after k iterations. By taking 






one easily calc, dates that the the Ith Fourier coefficient of error at iteration step k is given 
explicitly by 

(Note that the normalization of u imphes « = 0.) Taking <- = p/Ax < 1, we have that 

for each 1 = 1,2 J - 1. Moreover, the high frequency components of the error are 

contained in the Fourier coefficients ef with I » J/2. From above, we see that these high 
frequency components decay at a rate on the order of 

|1 - 2<r\. 

Therefore, these coefficients decay at a rate that is independent of the number of grid 
points J. On the other hand, the low frequency components of the error are contained in 
^th l « 1 or with l « J — 1. These components decay at a rate on the order of 

1 - ^T^C 1 " *)• 

Therefore for this example, doubling the number of grid points J can necessitate many 
mtre thaii twice the number of artificial time iterations to achieve a given error reduction. 

The basic idea of multigrid is the Mowing: One attempts to apteel 
components of u on a coarse mesh. This can be accomplished by a 

coarse grid information is then included into the current guess on the fine gnd, and l the 
hieh frequency components of the error can be rapidly damped there by using 

the one above. Normally, the divide and conquer strategy is employed 

in multigrid. That is, a nested sequence of coarse meshes is utilized to rapi y cap 
frequencies of the solution. 


To illustrate the multigrid approach, we will apply an ideal two grid multigrid strategy 
to example (2.1) above. Given the current approximation u* to u, we have that 


Dju + f = Dj(u k + c fe ) + f = 0 


or, 

(2.3) Dje k + r* = 0, 

where the residual is defined by r* = Dju k + f, and Dj is the J x J difference matrix 
defined above. The solution to the residual equation (2.3) is approximated by solvin g 

D J(2* k J/2 + r J/2 = 0 

where, 

r J/2 = 

I(J-+J/ 2 ) denotes an injection operator that takes vectors from R J into the lower dimen- 
sional space R* 7 / 2 . (Solving the coarse grid equation above exactly is what defines an ideal 
two grid strategy. In practice, multigrid is nested. That is, a nested sequence of lower 
dimensional multigrid iterations is applied to a nested sequence of coarser grid residual 
equations.) The fine grid error is then approximated by 

~ J(J/2-J)«J/2 = e cj» 

where denotes an interpolation operator that takes vectors from R ' 7 / 2 back into 

R J . e k g is called the coarse grid correction. One expects that e k g will agree well with the 
true error in the low frequencies. However, it in general will not agree well with the true 
error in the high frequency regime. To capture the high frequencies, we may for example 
use the artificial time scheme above, e.g., for n = 1, 2, . . . , v 

(2.4) e n = e n_1 - p(Dje n - x + r fc ), 


with 



(This step of multigrid is often referred to as the smoother. The term smoother is in fact 
a misnomer. In actuality, this step is used to capture the high frequency components of 
the error.) Finally, since 


u = u* + e* « u* + e*', 


we update u* via 


u fc+1 = u* + e*'. 


This defines one cycle of a simple two level multigrid technique. We should remark that 
the update step above is equivalent to setting 

u* +1 = SJ(u‘+4,f), 

where Sj(v, f) denotes v iterations of the basic iteration scheme (2.2) on grid J. 
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§2.2 Ideal two grid analysis for the example first order one sided scheme. A 
multigrid strategy requires three basic ingredients. They are: 

(i) The so-called smoother; see (2.4) above. 

(ii) The fine to coarse grid residual injection operator; I(j_ j/2) above. 

(iii) The coarse to fine grid error interpolation operator; /(J/ 2 -J) above. 

The performance of a multigrid algorithm relies heavily on the choice of these ingredients. 
In this subsection, we analyse the convergence rate of the strategy outlined above taking 
the artificial time relaxation technique (2.2) as the smoother, and 


(2.5) 


(■f(j— » j/2) r )i = •O’V + r 2i— »)i 


as the fine to coarse grid injection operator, and 

(2.6) (J(J/2-*J)*)2j = € i 

(f(J/2— j) c )2j-l = £ J> 


as the coarse to fine grid interpolation operator. In the context of finite volume, these inter- 
grid operators seem most natural. A multigrid iteration step above is written symbolically 

as 


(2.7) ->/3) r *) 

u‘+‘ = SJ(u‘ -Mo- 
using the fact that u = Sj(u, f), we may write (2.7) in terms of the error e k as 


( 2 . 8 ) 




= f(J/2— J)Sj/2(°> ) 

e fc+1 = S v j(e k - e* ,0), 


where again we normalize u such that — u j) — Sj=i € j ®* 

Unfort imately, unlike the single grid example above, the two grid multigrid algorithm 

(2.8) does not completely decouple into Fourier modes. Following (3], a J-penodic grid 
function v can be decomposed into a Fourier series on the odd grid points, and a Founer 
series on the even grid points. The explicit formulae are 


(2.9) 



J/2 

5 > e "^ ,(i+1)/2 

i= 1 


1=1 


if j = l,3, 


if j = 2,4, 


,J-1 
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In fact, if vj = 2ZiLi ^l e 2 - r a simple calculation will reveal that 


( 2 . 10 ) 


«l = -y (®l “ ^+J/2) c 2 ^', 
i>l = -y («l + vj+J/2) , 


for each 2 = 1,2, ... , J/2. 

(2.9) and (2.10) allow us to examine the degree to which the coarse grid correction 
€ cg approximates the true error e k . Writing 




or equivalently 


V.%. 

II 

-A. 

1 J / 2 

J-i 

we compute that 


-if.***). 



-*)■ 

Therefore, on the coarse grid, with j running from 1 to J/2, we have 


= 5 ((Bj £ ‘) 2 ,-, + (A,.* 

)*/) 


= L-fv* (1 «W)Sf«- 

v /J/2^2Az v ,l 

Solving the coarse grid problem exactly leads to 

S 

j/a(°» -^(J—j/2)-Dj«*) = 

i J / 2 , 

1 


and we find, after applying the coarse to fine grid interpolation (2.6), that the coarse grid 
correction is given by 

J/2 


<« 4 fc — ( Afc -^5 £*•-*** 
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Making us, of the formulae in (2.10), we can write the coarse grid correction in terms of 

{8f)i Li 


(2.11a) 

(2.116) 


(i^)l = \( 1 + e + 

(Z = 1,2, , *7/2). 


Notice that for l * 1, the coarse grid correction coefficients given in (2.11a) agree well 
with^he* true error coefficient, modulo a high frequency mmponent. The 
follow, from (2.11b) with 1 + J/2 « J. That is, modulo hrgh frequency P^rion, ‘ he 
coarse grid correction accurately predicts the low frequency components of the error. 

Finally after applying v iterations of the smoother (e* + -5j(( — e efl ,0)), one cycle 
of thfs idTri mudtigrid Nation yields an error which satisfies the coupled system 


/ i* +1 \ l( <(1 - e 3-‘) 

< 2 - 12 > UsJ-sU 


-o|'(l + e- 2 5“) 


°iW 1 + e 


■* 5 “) \ f /! 

-*•)) U+. 


)■ 

l+J/2 / 


where 


a\ = ((1 - <r) + tre*? 1 )*. 

We calculate the eigenvalues of the amplification matrix above. 

Aj = 0 

\ 2 = ]-(a " (1 - e -2 ^ 1 ) + ar +J/2 (l + e ***))• 

In the special «se when ^ = 1 and <r = 1/2, we have that A, = 0. Therefore tM»™ ltigri d 
sdiemerields the enact solution u after at most two cycle, regardless of it. uutial gums. 
Mo« generally for any v > 1 and any fixed 0 < <r < 1, one can show, mdependently of 
number of grid points J, that |A a | is strictly 1... than one. So in the generd lhe 
of convergence of this ideal multigrid scheme does not degrade with grid sire. 

A. mentioned above, the ideal scheme is not practical. One doe, not wish m sdv. Ae 
come grid problem exactly. Nevertheless, the analysis above does indicate that the fully 
nested multigrid approach should converge rapidly in a grid sise independent manner. 
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§3.1 A hybrid of multigrid. Consider for the moment the possibly nonlinear differential 
equation 


Q 

+ /(*) = °» 

together with given boundary conditions. Assuming that this equation has a solution, 
Newton’s method is developed by setting u = u k + e fc , where t t k is the current guess to the 
solution, and expanding g(u k + e*) = g(u k ) + g'(u k )e k H . Therefore, 

0 = ^ 9 ^ + ^ + + /(*), 


dx 


and we let e* m solve the linearized equation 


^(j’(«*)‘L) + ^s( «‘) + /(*) = o, 


to update 


u fc+1 = u k + e*. 

nm 


This idea extends equally well to finite diiference schemes associated to this differential 
equation. That is, if we seek the solution to the finite difference scheme 


Dj(g] u) + f = 0, 

where Dj(g\ u) denotes a finite difference operator consistent to d x g(u), we solve 

u*))4m + Dj(g; u k ) + f = 0, 

where d u {Dj{g-, u 4 )) denotes the Jacobian matrix of Dj(g\ u) at u = u*, and update 

i 

u T c nm* 


Newton’s method has two main drawbacks: (i) Each iteration requires the inversion 
of a large linear system, (ii) The Jacobian matrix of Dj(jj u)) may be quite complicated 
to compute analytically. In fact, for many modern finite difference schemes, the function 
u)) may lack the necessary smoothness to gain the quadratic convergence offered by 
Newton’s method. We propose the following simplification: Suppose we wish to solve 

Dj k% \g ; u) + f = o, 

where Dj *\g", u) denotes & (possibly high order) finite difference operator consistent with 

3»(<7( u )). ^ ^S ,O V(t0)« denote a (possibly low order) finite difference operator that is 
consistent to 5*(p # (u)e). Rather than solving the correct linearization of the fini te difference 
scheme above, we solve (or approximate the solution of) 

(3.1) D$'“V(u‘))ii_ + Bf'te u‘) + f = 0, 
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to update u fc+1 = u fc + e*~. Certainly there is no reason to expect quadratic convergence 
from this approach. Themfore, rather than solving (3-D .«<*. - W^dmate e~ by 
applying multigrid to (3.1). Specifically, suppose that u‘ is known. Then 

(3.2) 

Step 1: Compute i£(u fc ) = Dj \g~, u ) + f* 

Step 2: Apply n > 1 cycles of multigrid with initial guess 0 to approximate 

the solution to 

bS , 'V(»‘))‘~+ b ( u ‘) = 0 ' 

calling the result 

Step 3: Update u *"*” 1 = u* + c m p* 

Step 4: (Optional) Steps 1-3 may not capture the high frequency components 

of u well. Therefore, applying a few iteration of a smoother consistent 

to X>J f) (s; u) + f = o may be advised; (see below). 


83.2 Convergence of the hybrid multigrid approach for a model Problem. Wbetber 
iterates coming from (3.2) converge or not depends in a cruel way on he choice of both 
D (,o) and M W) . We again only perform ideal two grid multignd analysis to the linear 
problem g( u) = u. The problem is made somewhat more interesting by taking an upwind 

third order scheme for D j ^ 

(3.3a) (DS‘°“)f = - “!->) " 6^ ( “ i+1 _ ^ + 3 “ i_ ' " 

and a first order upwind scheme for Dj * 

(3.36) 




together with periodic boundary conditions. (Not. that (3.3.) is third order in the sense 
of cell averages.) We still require that /( z i) ” 0 for solvabait y* ^ 

Before studying (3.1) & (3.3) hybridized with multigrid, observe that if e k = u - u 
(where again u is normalized such that = 0), then e~ (see 3.1) satisfies 


(3.4) 


x nm' 


where the amplification factor V*i P ven 

(— i sin (flt) — l/3(cos(#t) — 1)(1 ~ e *)) 
= (t-e«0 


= 


«. = 7 '- 
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The error reduction r&te expected from algorithm (3.1) (or equivalently algorithm (3.2) 
with Step 2 solving exactly) is on the order of maxi<j< j_i |1 — if>i\ ~ 0.53 and the ma ximum 
is attained in the mid frequencies. For the high frequencies l fa J f 2 we have that |1 — ij>i\ fa 
1/3, and for the low frequencies IfaloilfaJ — lwe have that |1 — -0i| « 0. So we see 
that the low frequency components of the error are captured quite rapidly whereas the 
mid to high frequencies decay less rapidly. Step 4 can therefore be utilized to knock down 
the mid to high frequency components of the error without affecting the low frequency 
convergence. 


The results of Section 2 make the analysis of the hybrid multigrid algorithm (3.2) 
applied to the difference operators in (3.3) an easy matter. Since — e* is the error of 
Step 2, using (2.12) we find that 



So if we denote the 2x2 multigrid amplification matrix above by M \ a ' v,v ^ , and define 

- (+* 0 ^ 

' \0 i> l+J/2 J ’ 

then referring back to (3.4) we find that the 2x2 matrices 


(3.5) Ai= (/-(/- M} <, ™ ) )h) , (Z = 1,2,... , J/2), 

define the amplification matrix to the full (Steps 1-3) hybrid algorithm (3.2). The spectral 
radius of the hybrid amplification factor is plotted as a function of a in Figure 1. (Recall 
that <r is the ratio of the artificial time step size p to the grid size Ax in the approximate 
linearization multigrid smoother.) Note that for a = 0.40 we actually achieve better error 
reduction using multigrid in Step 2 than if we solved (3.1) exactly. 



Figure 1. The spectral radius of diag(Ai) plotted against a. 
\y - 2, rj = 2.) 
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Notice the sensitivity of the error reductxon rate to * m Figure 1. This ? 

cates the relationship of multigrid error reduction to linearization wave speeds. T ’ 

for nonlinear problems where the wave speeds are not constant, pn^itive 

systems with multiple wave speeds), we expect a poor convergence rate from p 

algorithm outhned above. This problem will be addressed in the next section. 

84. Application to the one dimensional Euler equations. The partial different^ 
equations that govern the flow of a compressible and invisad gas m a quasi one dimensional 

expanding duct are 

J^) + = 0 


dx' 

l(puA) + ^ ((pn 2 + P)A) - p£-A = 0 
| j{peA) + ^ ((pe + p)uA) = 0, 


where p is the fluid’s density, tx its velocity and e its total energy per unit mass. The given 
function A = A(x) defines the cross sectional area of the duct as a function of positio 
along its length, p represents the fluid’s pressure and is given by the equation of state 
/ jw e _ u 2/2). We take 7 = 1.4 in the results presented below. We seek a 
steady state solution to problem on the interval 0 < x < 10, taking a supersomc inflow 
boundary condition at x = 0: (p, u,p) = (0.502,1.299,0 3809); and a subsomc outflow 
condition at x = 10: (p,u,p) = (0.776,0.5003,0.7475). The cross sectional area given 

by i4(*) = 1.398 + 0.347 tanh(0.8x — 4). 

To demonstrate a need for a technique that speeds convergence to steady state, we seek 
a steady state to this duct example by applying the third order accurate 
scheme from the introduction together with artificial time relaxation. Specifically, 

iter&te k = 1,2, . . . 

q‘ = q‘-> - pK* 1 ' V). 
where the residual is given by 

(«< M >( q))i = (h*(q,+i/a;qj + i/j.qj+i/z) ~ h s(<lj-i/ri q i-i/i' < 'j-i/i)) + 

The compression factor in the slope limiting above as well as in ^ other examples 
is taken to be 3. We initialize the iteration by a linear interpolation of ‘J' 1 " 1 " ! 

in the conserved variables. The relaxation factor p is taken so th. the CFL munW(th. 
ratio of time step size times the fastest wave speed to the space step size) 0. 
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Figure 2. §log(^j (q k ))*) plotted against iter*- 

tions for the artificial time scheme. 

use 64 grid points on the interval [0,10]. In Figure 2, we plot the base 10 log of the L 2 
norm of the (first component of the) residual versus the number of iterations. 

To apply the multigrid algorithm (3.2), we need to decide on the approximate lin- 
earization of the third order scheme used above. We take for -D (io *(ff'(q))e the following 
first order scheme 


Ax c i> c i+i) ” ^^(qj-i^; €j_i, Cj)), 

where 

h fiq8^ q j +1 /2; »«/+*) 

= 2 + € j) — l^qg(<l;+l/ 2 )|(Cj+l — «j)). 

As usual, |dqg(<i)| denotes the matrix -R(q)| A(q)|£(q). The obvious choice for a smoothing 
iteration scheme is 

4 " = «?•«*«) 

- •O’fc-'/* 4-«4)) + M«u>4 + (^“’(q)),), 

where p is fixed according to the largest wave speed of the variable coefficient problem. 
However, from our analysis in Section 3, we can expect this approach to yield very poor 
convergence rates. We base this conclusion on the following observation: The choice of p 
for the smoother above must be based on the largest eigenvalue of 0 q g. This is needed for 
stability. Since the flow contains a supersonic region, the eigenvalue u + c (c is the local 
speed of sound and u > c is assumed) is quite large. Therefore, p must be taken quite 
small. On the other hand, near a sonic region, u - c is quite small. Therefore with this 
simple smoother, the hybrid multigrid performance in sonic regions can be predicted by 
the results in Figure 1 with <r near zero. There, the error reduction rate is near uni ty. 
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For this reason, we propose characteristic time stepping. This is accomplished by a simple 
modification of the smoothing iteration above: 

e* +1 = tj - pR(qj)|l/A(<lj)l i, (qj)^^0 1 « < rf(qi+ 1 / 2;e i ,€ >+ 1 ) 

- + M«u)4 + (^“’(q));). 

where now, p is taken to be p = 0.4 Ax. 

We implement algorithm (3.2) to this duct test problem, again using 64 grid points 
on the interval [0,10]. Multigrid Step 2 is done using every grid level 2 1 , 2 2 , 2 3 , 2<, 2 5 , 
2®, with v = 4 and rj = 1 . Step 3 is relaxed somewhat, specifically q fe+1 = q +5 e mp , 
because of the presence of the supersonic-subsonic shock. Around the shock, the slope 
limiting algorithm introduces a small amount of numerical viscosity and this viscosity can 
introduce an instability into the approximate Newton iteration; (see equation (3.4)). Step 
4 of the hybrid algorithm is also utilized with 2 characteristic time stepping iterations 
applied to the third order residual. 


0.0000 
-1.0000 
-2.0000 
-3.0000 
-4.0000 
-5.0000 

a 0000 

0.0000 5.0000 10.000 15.000 20.000 25.000 

Figure 3. Log residual versus hybrid multigrid cycles, again ap- 
plied to the duct problem with 64 grid points. 

Comparing Figure 2 with Figure 3, we see that the multigrid algorithm reduces the residual 
of this test problem in 100 times fewer iterations than the artificial time scheme, whereas 
each multigrid cycle costs on the order of only 4 times the work. Figure 5 depicts the resid- 
ual reduction from the hybrid multigrid scheme on 5 grids; from level=4 (16 grid points) 
to level=8 (256 grid points). This last figure demonstrates the virtual grid independence 
of the hybrid algorithm. 
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Figure 4. The computed density for the duct test problem using 
64 grid points. 


# levels ( 4 5 6 7 fl ) ~ ( +.x od x ) 



0.0000 5.0000 10.000 15.000 20.000 25.000 30.000 35.000 40.000 45.000 

Figure 5. The log residual versus hybrid multigrid cycles for the 
duct test problem using 16 to 256 grid points. 
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